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1. Introduction 

Delamination in the form of cracking or separation between plies in an advanced fiber 
composite laminate has been a problem of major concern. Delaminations are recognized to result 
from manufacturing/processing defects as well as from high, local interlaminar stresses at 
geometric boundaries or under impact. The presence of a delamination and its growth may have 
significant effects on mechanical response and load-bearing capacity of the composite laminate. A 
delamination may reduce local stiffness structure and, in many cases, lead to premature failure of a 
composite structure due to local buckling and/or appreciable crack growth. This situation is most 
critical when a delamination in a composite laminate is subjected to out-of-plane loading, such as 
bending and torsion. A rigorous study of the delamination problem for a fiber composite under 
out-of-plane bending is recognized to be difficult. The difficulty is caused by various complexities 
involved , such as inherent multiphase material system, significant anistropy of the composite, 
discontinuous mechanical properties through the laminate thickness direction, a large number of 
lamination parameters, and associated crack-tip stress singularity. 

Extensive research efforts have been made on the delamination mechanics problem of a 
composite laminate subjected to in-plane tensile or compressive loading, for example [1-51. 
Interlaminar stresses are found to be dominant in governing the initiation and growth of the 
composite laminate under an in-plane nominal stress. It is recognized that the interlaminar stresses 
are more prevalent in a fiber composite laminate subjected to out-of-plane bending [6] than that to 
in-plane stretching. Consequently, a delamination in a fiber composite laminate subject to bending 
would experience more significant influence by the presence of the interlaminar stresses. Accurate 
determination of the magnitude and distribution of interlaminar stresses becomes critical. Attempts 
on analytical evaluation of the delamination behavior in a composite laminate have been extremely 
limited [ 7 ], because of the difficulties in accurately determining the complex interlaminar stresses 
and deformation. It is not surprising to recognize this difficulty since interlaminar stress and 
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deformation are generally three-dimensional in nature. The local crack-tip response associated with 
a delamination requires simultaneous evaluation of the coupled mode-I, II and III fracture 
mechanics parameters. However, in order to have a clear understanding of the delamination 
problems in a composite material under bending, an accurate solution obtained based on the 3-D 
(or quasi-3-D) laminate elasticity theory and the recently developed composite fracture mechanics 
analysis is essential. Both quasi three-dimensional solutions and fully 3-D solutions are 
considered to be critical to the advancement of analysis and design of composite materials in the 
future. Also, with the aid of the recently developed laminate elasticity for fiber composites, it 
becomes possible to look rigorously into the delamination problem in a composite laminate under 
bending. 

The objectives of this paper are to: (1) develop an asymptotic solution for a composite 
laminate subjected to out-of-plane bending, (2) construct advanced singular finite elements in 
conjunction with the developmenmt of nonsingular elements for this bending problem, (3) evaluate 
the delamination failure mechanics parameters and the subsequent modes of fracture. Both 
advanced analytical methods and advanced computational analyses are conducted to realize these 
goals. A parametric study was also conducted to evaluate the influences of various lamination 
parameters on the delaminated composites. 

In the next section, basic formulation of laminate elasticity is given and governing partial 
differential equations of the delamination problem are established. Also, the general solution 
structure is determined by an eigenfunction expansion method. In Section 3, advanced numerical 
methods are introduced to formulate singular, local crack-tip elements and the nonsingular elements 
in the surrounding area in a delaminated composite. The fracture mechanics parameters, such as 
delamination crack-tip stress intensity factors and energy release rates, for the problem is addressed 
in Section 4. Numeral results are obtained for various cases of interest in Section 5. Important 
problem parameters, such as the size of the singular element, the effect of delamination length and 
fiber orientation as well as the laminate layup, are studied to examine the fundamental nature of the 
solution method and the delamination crack-tip field. The effect of delamination location, i.e.. 



3 


edge delamination and center delamination, on the local fracture is examined in Section 5. 
Conclusions derived from this study is given in Section 6. 

2.1 Basic Formulation 

Consider a composite laminate composed of fiber-reinforced composite plies subjected to 
an applied nominal end bending moment per unit width, M^, as shown in Fig. 1 . Let V denote the 

interior of a composite laminate of length 2L and a cross section B with the lateral boundary 9B. 
The Cartesian coordinates are chosen such that the z-axis is the generator of B. The end cross 
sections are denoted by and respectively, and located at z = -L and z = L. The constitutive 
equations of each individual composite ply are denoted by generalized Hooke's law in contracted 
notation as 


£i = Sjj Gj, (i, j =1,2,3,.. .,6) 


( 2 . 1 ) 


where repeated subscripts indicate summation, and Sjj is the compliance tensor. The engineering 

9w 

9u 9v 


strains, ej, in Eq. (2.1) are defined by 

9u 9v 

e l = £ x - e 2 - e y - 

dx dy 


e 3 = e z = 


9z 


9w 9v „ 9w 9u 

e 4 = 27y Z = - + i? e 5 = 2Yxz = ^. *6 = 2Yxy = ^ 


( 2 . 2 ) 


where u, v and w are displacement components. The stresses, cj, are defined in an analogous 


manner. 

The composite laminate is assumed to be sufficiently long that in the region far from the 
end, the end effect is neglected by virtue of Saint- Venant's Principle. Thus in this generalized 
plane deformation problem [ 8 ], stresses in the laminate are independent of the z-coordinate. 
Under these assumptions, equations of equilibrium without body force lead to 
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3q x ^T X y _ 

3x 3y 
3t XV 

— M + — ^ = 0, 
3x 3y 


d^xz + ^ T yz _ 

3x 3y 


Integrating Eq. (2.1) with the aid of Eq. (2.2), one can obtain u, v and w as 

z 2 3D 3w 0 

u = - y — + (S 5j aj - z + U 0 (x,y) 

z 2 3D 3w 0 

v = - y — + (S 4 jaj - -^) z + V 0 (x,y) 

w = Dz + W 0 (x,y) 

where D = S 3 j O) and U D , V 0 , and W 0 are arbitrary functions of x and y only. 


(2.3) 


(2.4) 


(2.5) 


Following the procedure in [ 9 ], it can be easily shown that general expressions for the 
displacements and the stress component a z have the following form: 
u(x,y,z) = y Ai S 33 z z -A 4 yz + U(x,y) 

v(x,y,z) = y A 2 S 33 z^ + A 4 xz + V(x,y) (2.6) 


w(x,y,z) = (Ajx +A 2 y + A 3 )S 33 z + W(x,y) 
a z = Ajx + A 2 y + A 3 -S 3 jaj/S 33 , (j=l,2,4,5,6) 

where A] and A 2 pertain to the external bending on the composite laminate, acting in the x-z and y- 
z planes. The A 3 characterizes uniform axial extension of the composite, and A 4 , the relative angle 
of rotation about the z-axis. 

The unknown functions, U, V, and W, depend on x and y only, and can be easily shown 

to obey the following relationships: 
d\J ~ 

— - = Sij Gj + S 13 (Aix + A 2 y + A 3 ) 

3x J J 

3V ~ 

— = S 2 j cjj + S 23 (A ix + A 2 y + A 3 ) 

3W ~ 

— = S 3 j aj + S 53 (Aix + A 2 y + A 3 ) + A 4 y 


(2.7) 
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aw 

9 y 


= §4j Oj + S43 (A|x + A2y + A3) - A4X 


au av ~ _ /A . . , 

— + — = s 2j <*j + S63 ( A l x + A 2y + a 3>’ 

ay ax j j 

(j = 1,2, 4, 5 and 6) 


where Sjj = Sjj . Sj3 Sj3/S33 t ( 2 . 8 ) 

Applying the interface continuity conditions, = u^ m+ ^’ between the m-th and (m + 

l)th plies leads to the following relationships for Aj ; 


A i (m) S 3 3 (m) = A i ( m+1) S 3 3 (m+1) > (2 - 9 ^ 

A 4 (m) = A 4 (m+1) - (i = U. 3 ). 

Traction-free boundary conditions on the lateral boundary aB satisfy the following 
expressions: 

CT x n x + X xy n y = 0, 

^xy n x + °y n y = 0, 
x \z n x + T yz n y = 0, 


( 2 . 10 ) 


where n x and ny are outward unit normal of aB. 

The end loading conditions on the laminate cross-section under the present bending M x , 

have the form obtainable from statically equivalent loads as 


^xz dxdy = 0, ^yz dxdy = 0 

^z dxdy = 0 
T xy dxdy = ^x 

X yz dxdy = 0 (i = 1,2) 

% ( X yz x - X xzy) dxdy = 0 


( 2 . 11 ) 
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2.2 Governing Partial Differential Equations 

Introducing stress functions F(x,y) and y (x,y) such that 
0 2 F a 2 F a 2 F 


a x = 


Sy r °y = ax 2' ‘*y = - 3x3y ' 


tyv = 


X XZ = 


By 

By’ 


B\\r 


Tyz = - Bx 


a i2) 


the equations of equilibrium can be satisfied identically. Eliminating U and V from Eqs. (2.7) by 
differentiation, one obtains the following system of coupled governing partial differential equations 
for this problem: 

L 3 F + L 2 Y = - 2 A 4 + A^34 - A2^35, 

L 4 F + L 3 \|/ = 0. (2.13) 

where L 2 , L 3 , and L 4 are linear differential operators defined in [ 9 ] 

~ B 2 ~ B 2 ~ B 2 

L 2 = S44 —z - 2S 45 r - t - + S55 


0 x 2 


0 x 0 y 


0 y 2 


~ 03 ~ ~ 03 ~ 03 ~ 53 

l 3 = S 24 ^3 + (S25 + S 46) ^ - <S 14 + S 56) — + S 15 ^-3 


(2-14) 


~ 04 ~ 34 ~ ~ 04 04 ~ 04 

l 4= S 22 ^ ■ 2S 26 ^- + (2S 12+ S 46) ^ ' 2S 16 ^ + S 11 ^ 

2.3 General Solution Structure 

The governing equations Eqs. (2.13), are coupled linear partial differential equations with 

constant coefficients related to elastic constants of each composite lamina. The complete solution 

for Eqs. (2.13) consists of a homogeneous solution and a particular solution. Lehknitskii has 

shown [ 8 ] that the homogeneous solution for the above system of governing partial differential 

equations has the general form 

6 

F(x,y) =X F k ( x + Hky)> 
k=l 


(2.15a) 
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V(x.y) =y\kFk(x +Hky) (2.15b) 

k=l 

where the prime (') in Eqs. (2.15) denotes differentiation of the function Fj c (x +(i| c y) with respect 
to its argument. The coefficients are the roots of following characteristic algebraic equation: 

400 /2OO - /3 2 (M-) = 0 


where 

*200 = S 55 |i 2 . §45 p 2 + §44 (2.1 6 ) 

400 = S 1 5 ^t 3 - (S 14 + S56) it 2 + (S25 + §46)11 - §24 
400 = §llli 4 - 2 § 16 l l3 + ( 2 Sl2 + S66)ii 2 - 2 S 26 H + §22 

It can be shown that cannot be real and have to appear as complex conjugates. Thus are 
analytic functions of the complex variables = x +|i} c y. The general solutions for stress and 
displacement are constructed in the form as 

°i = a i (h) + a i (p) , (i= 1,2,.. ..,6) 

(2.17) 

uj = uj (h) + uj (p) 0 = 1.23) 

The homogeneous components of the stress and the displacement solutions can be expressed in 
terms of Fi c (Zj c ) as 


a v (h) = 


Y *W Z k>. a y O) - £ F' k (Z k ) 
k=l k=1 


yz^) ^k^k^k^’ ^xz^ ^ ^ ^k^k^k^k^ 

k =l 1 tl 

6 


T X y (h) - 


y ^ F k< z k) 

k=l 


(2.18) 
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6 6 

„(h) = Y p k F k (Z k ), v(h) = Y, IkFkW. 

k=l k=l 

6 

w(h) * X 'k F k< Z k>- 

k=l 


(2.19) 


where pj^, and are known constants with their expressions defined in [ 8 ]. 

The aj(P) and uj(P) are particular solutions for the problem, depending upon the loading 

mode, composite lamination parameters and geometric variables of the composite. The particular 
solutions for aj and Uj can be found by various analytical and numerical methods, dependent upon 

the complexity of the problem. 

An appropriate form of Fj ( (Zj c ) has been suggested in [ 9 ]. 

F k (Z k ) = c k Z# + 2 ) / [(8 + 1) (8 + 2)] (2.20) 

where C k and 8 are complex constants to be determined later. Substituting Eq. (2.20) into Eqs. 

(2.18) and (2.19), one can express oj (h) and uj 0 1 ) in terms of 8, C k and Z k . 

The homogeneous solutions for the stress and displacement shown in Eqs. (2.18) and 

(2.19) are required to satisfy the near-field boundary conditions and interface continuity 
conditions, leading to a standard eigenvalue problem in the form of a transcendental equation for 
determining the values of 8. After the eigenvalues, 8, are determined, the stress and 
displacement may be expressed as 


°i ^ Z ^ d m ^im(Zk’^m) 
m=1 k=l 

oo °o 

uj (h) - Y V d m gjm(Zk^m) 
m=1 k=l 


( 2 . 21 ) 


where the d m are unknown constants and fj m and gj m denote the known eigenfunctions 
corresponding to the m-th eigenvalue 8 m . 
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The particular solutions for the stress and displacement in the present delamination problem 
of a composite laminate under out-of-plane bending are extremely complex and difficult to be 
obtained analytically in a closed form due to complex geometry of the crack, the laminate 
parameters and the out-of-plane loading mode. An advanced numerical method is developed in 
conjunction with the asymptotic solution obtained here to give the full-field solution for the 
problem. 
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3. Numerical Methods 

Owing to the complex nature of the problem, a closed-form full-field solution for the 
delamination problem in a fiber composite laminate under out-of-plane bending is difficult to 
obtain, especially the determination of the particular solution. Advanced numerical methods are 
introduced here to account for the various difficulties and to obtain the full-field solution. Based 
on the asymptotic solution determined in the previous section, a singular delamination crack-tip 
element is formulated. Also, based on the laminate elasticity field equations, a compatible non- 
singular surrounding element is constructed for the bending problem of a delaminated composite 
laminate. 

3 . 1 Near-Field Element Formulation 

In the previous studies [10,1 1,12], the detailed nature of the composite delamination crack- 
tip singularity has been examined. For a delamination with a closed crack tip with frictionless 
crack surface contact, an inverse square-root singularity has been determined [12]. Also, for an 
open delamination crack tip and for a crack with a small amount of crack closure, asymptotic 
solutions with an inverse- square root stress singularity are also approximated [11]. Based on these 
results, the delamination crack-field solution structure can be properly constructed. In order to 
include the local stress singularity, the element formulation introduced by Stem[13] is employed. 
The element is a six-node triangular element in which the assumed displacement field contains 
terms proportional to the square root of distance along lines emanating from the crack tip. The 
element is fully compatible with the neighboring isoparametric quadratic elements. 

Nodal arrangements and element configurations for the present problem are shown in Fig. 
2 with the node 1 being singular. Transformation is needed for the triangular polar coordinates 
(p£) to the global Cartesian coordinates (x,y,z). Coordinates of any given point in the element are 
given by 


x = xi + p { ( X 2 - xi) + £, ( X 3 - X 2 ) }. 


(3.1) 
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where x^ (i=l,2,3) are global Cartesian coordinates of the vertices of the singular element. 

Furthermore, the relationships between the spatial polar coordinates (r,0) with the origin at node 1 
and the triangular polar coordinates (p,5) are 

r = R (5) p, 6 = 0(5). (3-2) 

where R (5) = { (x 2 - x^ 2 + (y 2 . y^ 2 + 2 5 [(x 2 . (x 3 . x 2 ) + (y 2 . y^ (y 3 . y 2 )] 

c 2 r , , 2 , , .2, ,1/2 
+ 5 K x ! ■ x 2 ) + (y 3 -y 2 ) 1) 


0 (0) = tan 


■ i (y2-yi> + Uy 3 -y 2 ) 

(x2 - Xj) + 5 ( x 3 - x 2 ) 


Interpolation functions Nj for the displacements u, within the element are expressed as 


2-2 (1_X,) p X 

N = 1 + P ' K 

1 2 (1-X) . j 2 (1-X).i 


) 


N = 


2 (1-X) 

2 P(1 - y ^xT 

>(1-X) 


- pX [25(1-5) + 


H 


-l 


,d-X) 


-i 


n 4 

N 5 = 4pX5 (1-5) 

N, = -^- + 2p^-S— 

6 2( 1 ’X)-1 2(!-X)-l 

The displacements u^ = {u,v,w} can be decomposed into the following: 

6 

u = { pX Nj (1) (5)+ n. (2) (p,5)} u. 

i=l 

Along 0 = 0 O , i.e., the angle of the crack surface, the displacements are approximated by 

6 

i=l 


(3.3) 


(3.4) 


(3.5) 



12 


where 


. ,i ( x 2 • x i) tan0 o - (y 2 - yp 

S|n_n — 

0 (yg - y 2 ) - ( x 3 - x 2 ) tan0 o 


Thus the singular terms in the stress expression are 

6 


G = r 


X-i 


Qi(^o) Hi 
i=l 


(3.6) 


where Qj is a matrix with elements related to the derivatives of shape functions and elastic 
constants. 


3.2 Far-Field Element Formulation 

elements in the solid composite modeling are formulated based on generalized 
plane deformation theory and the minimum potential energy principle. Since no stress singularity 
is involved in the far field, formulation of the displacement-based element in a composite laminate 
under bending requires less effort than that given in the previous section for the construction of a 
singular element. A triangular element of the same nature has been proposed by Wang, et al[ 14] 
for a delamination problem in the composite laminate under in-plane stretching. To improve the 
rate of convergence of the solution, an eight-node isoparametric element is developed in this study 
in conjunction with the singular crack-tip element described in Sec. 3.1. A brief outline of 
derivations of the element stiffness and the associated loading vector of the far-field element is 
given here. 

The potential energy functional II p to be minimized for the current element formulation in 
the absence of body force is given by : 

= lJA m 2 C e dA - JJs T*^uds 

mz~ — a m ~ 


(3.7) 
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* 

where e is a strain vector, C is a 6x6 composite stiffness matrix, and T , the prescribed traction 

along s a of the mth element with an area A m . 
m 

For the present problem of a composite laminate under bending M~, the displacement field 
may be written as 


u = u o + U (x,y) 


(3.8) 


where 

T 

u = { u,v,w } , 

u 0 T= {-1/2 - A^yz, -1/2 + A^xz, (A^x-h^y+A^) S^z } 

U T = { U(x,y), V(x,y), W(x,y) } (3.9) 


in which U(x,y),V(x,y) and W(x,y) are unknown functions to be determined. 
Using an isoparametric element representation for U, V and W, one has 


U = N q (3.10) 

where q are nodal displacements and N is a matrix of nodal interpolation functions. The modified 

strain-displacement relationship for the composite bending problem has the form, 
e = e + B q 

~ ~o ~~ (3.11) 


where £ — {£x’ e yiYyz>Yxz’Yxy}> e 0 - (A^x-h^y+A^lS^g, A^x, -A^y, 0}, and B is the 

differential form of N. Since N is a function of x and y only, it is convenient to rewrite Eq. (3.1 1) 
in a simplified form as 

*5x1 = B 5x24 q 24xl (3 - 12) 


_T 

where £ 


= { £ x , £y, Yy Z , y xz , Yxy h and B is the reduced form of the aforementioned gradient 


matrix B. 



The stress-strain relationship for a composite lamina under pure bending can be written in 


modified contracted notation: 


a= C (e - £q) (3.13) 

where 

a ^ = ( <*x, a y, x yz, X xz, X xy } 

_ y 

e 0 = t e 0i ' = t S 13 A 2 y ’ S 23 A 2 y ’ ' A 4 X + S 34 A 2 y ’ 

V + S 35 A 2 y, S 36 A 2 y }, 

0>x5 = { S 5 X 5 ] * (i,j=l,2,4,5,6) 

and Sjj is the reduced compliance tensor of each lamina. The stress component, C7 Z =0'3, is 

obtained from 

Og = AjX + A 2 y + - S^j C5j/S^^ (3.14) 


In terms of e and £q > the potential energy functional, lip, at an element level can be written 


as follows: 


n P = WA m ie T CedA-lJ Am ^eT C e 0 dA + i J^eJc^dA 


Js a T* T uds + C 0 


(3.15) 


-’m 


where C 0 is a constant related to £ 0 . 

Following the standard variational procedure of minimum potential energy, one can obtain 
the element stiffness k and the consistent loading vector Q as 
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£ e o 


d A + 



T * 

N X T ds 
m~ ~ 


(3.17) 


for a far-field element in the modeling of the bending problem of a delaminated composite laminate. 


3 . 3 Solution Scheme and Algorithm 

Based on a displacement-based finite element method, it is convenient to apply extensional, 
bending and torsional deformations, i.e., A\, A 2 , A 3 , and A 4 ,in the formulation. It is well 

known that, in general, bending of a composite laminate is accompanied by extension and torsion. 
Based on the superposition principle of linear elasticity, total deformation under mechanical loading 
can be expressed as a linear combination of the extentional, bending, and torsional deformations. 
Written in the mathematical form, one has 


A 1 S 33 


D 11 D 12 D 1 3 D 1 4 


M°° 

y 


A 2 S 33 


D 12 D 22 D 23 D 24 


M°° 

Y 






A 

(3.18) 

A 1 S 33 


D 13 D 23 D 33 D 34 


z 


- A 4 - 


- D 14 D 24 D 34 D 44- 


M~ 



where = M t = 0 and F z = 0 in the present case. The symmetry of the D matrix is assured 
by Betti's reciprocal theorem. 

By inverting the above equation, the induced extensional, bending and torsional 


deformations under applied bending are obtainable. 
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4. Delamination Stress Intensity Factors and Energy Release Rates 

Fracture mechanics parameters, such as crack-tip stress intensity factors and energy release 

rates for a delamination, are evaluated from the asymptotic solutions for interlaminar stresses 
°2’ a 6 anc * °4^ or G y’ X xy’ anc * X yi) anc * displacements u- along the plane of the crack. For a 

finite-dimensional composite laminate containing a delamination, the asymptotic stress and 

displacement fields can be conveniently determined by the aforementioned singular finite-element 
method. Along the delamination plane \ , the near-field stresses and displacements are 

approximated by using Eq. (3.7) and its derivatives as 

a=rp 5 R(^ o ;5) q (4.1a) 

u M(y q + u o , (4.1b) 

where R(^ q ; 5) is a matrix of derivatives of the shape functions M and ply elastic constants. The p 
and £, o are related to the global coordinates by a simple transformation given in [14]. 

For a partially closed delamination or a delamination with a very small size of crack-tip 

closure for which the simplified model [11] with an inverse square-root stress singularity is used, 
the stress intensity factors and strain energy release rates can be evaluated easily. Taking a and u 

along the crack plane J; = 0 (i.e., <}> 0 = 0), the asymptotic interlaminar stress and displacement 

along the ply interface can be written in a form as 

G. ~ A k r ' 1/2 (i=2, k=l; i=6, k=2, and i=4, k=3), (4.2a) 

u- ~ B,r 1/2 + u • (j=l, k=2; j=2, k=l, and j=3, k=3), (4.2b) 

J K U J 

where A^ and are obtained from the corresponding components of R(^ q ; 5) and M(J; o ) in Eqs. 
(4.1(a)) and ((4.1(b)) by setting t, = 0 and 8 = -1/2. Thus, the delamination stress intensity 
factors Kj can be easily determined by 

K i = V^T A. (i=I,II,III). (4.3) 

The energy release rates can be determined in a manner similar to that for K. as: 

G = G i + G 2 + G^ 
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lim 1 f8a 

= 5a->° 25a J ° 


3 

{ I (A, 
k=l * 


r' 1/2 ) [B. 


(5a-r) 1/2 3)dr 


r i ? 3 /8a-r x W 

lim J_ j6a z A B (— ) dr 

= 5a— >0 25a 0 k=l k k V ' 


(4.4) 


The term u Q j in Eq. (4.2b) is not included, because it does not result in any contributions to G- and 
G after integration. Integration of the singular integral Eq. (4.4) can be carried out explicitly. The 


strain energy release rates have the form 
3 


G =Ji G i -jOW a 2 b 2+ a 3 b 3 ). 


(4.5) 


We remark here that each term in Eq. (4.5) corresponds to an individual G- and that for a 


delamination with a finite length of crack closure, the first term in Eq. (4.5) is identically zero, i.e., 
G x = 0 . 
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5. Results and Discussion 

The aforementioned analytical and numerical methods have been applied to a delamination 
problem in a fiber-reinforced composite laminate subjected to bending to illustrate solution 
convergence and computational efficiency. A schematic finite element discretization for a half of 
the composite laminate cross section is shown in Fig. 2. A six-node singular element with 
eighteen degrees of freedom (3 DOF per node) around the crack tip, shown in Fig. 3, is 
constructed in conjunction with the eight-node, quasi three-dimensional, nonsingular element in the 
surrounding region. Compatibility between the two kinds of elements is ensured in their 
formulations (Secs. 3.1 and 3.2). The non-singular finite elements are refined near the crack tip 
region to take care of the high stress gradient. The influence of the size of the singluar element on 

convergence and efficiency of the present approach is demonstrated in this section. 

For illustrative purposes, numerical results for a [45/-45] s graphite-epoxy composite with 

two edge delaminations under bending is considered. (The fiber orientation 6 is defined as the 
angle between the fibers and the z-axis.) The composite has a geometry of b = 2 in., h = 0.25 in. 
and a/b = 0.35. Due to the out-of-plane loading mode and the geometric symmetry, only one half 
of the cross section is needed. Assume that each lamina is made of graphite-epoxy with the 
following elastic properties: 

E l = 19.5 10 6 (psi), E T = E z = 1.48 10 6 (psi), 

G LT = G LZ = 0 8 1()6 (P si )’ G TZ= 0 49 48 1q6 (psi)> 

U LT = V LZ = V TZ = °‘ 49 ‘ 

where L, T, and Z refer to the principal material axes along fiber, transverse, and thickness 
directions, respectively. 

For symmetric angle ply composites under pure bending, symmetry conditions for 
displacements with respect to the y-z plane lead to the following relationships: 

U(x,y) = U(-x,y) 

V(x,y) = V(-x,y) 

W(x,y) = W(-x,y) 


(5.1) 
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5.1 Size Effect of the Singular Element 

In the analysis of a delaminated composite laminate, the singular elements are required to 
cover a finite region at the crack tip so that the deformation in the singular-stress domain can be 
accurately determined. It is recognized that the singular domain is small in general and that the 
shape functions consist of only dominant singular terms and the terms for the nonsingular field. 
Thus sizes of the crack-tip elements are of significant importance in the solution accuracy and 
efficiency of numerical results. In this study, a total of twelve equal-size elements around the crack 
tip, shown in Fig. 3, are used to cover the singular domain. The size effect of the singular element 
on solution convergence and efficiency may be characterized by a non-dimensional parameter rj, 
defined as 

Aa linear dimension of the singular element 

^ - h lamina thickness 

In Fig. 4, the normalized stress intensity factors Kj (i = I, II, III) in the delaminated 
composite change as the size parameter rj reduces from 0.01 to 0.001. In all of the cases, the 
variations are less than 3 percent in the acceptable range. It indicates that the present solutions are 
converged and stable. Thus, in all of the later analysis, T| = 0.004 is used. 

5.2 Effect of Number of Elements around the Crack Tip 

In the singular element formulation, shape functions for the dominant singular, nonsingular 
terms (eq. 3.4) are quadratic and linear in 0, respectively. The minimum number of singular 
element enclosing the tip region, which would provide sufficient angular variations in the near 
field, is an important computational parameter to be determined for the solution convergence and 
efficiency. The 2jt/n is the size parameter for the solution convergency and accuracy. This effect 

may be clearified by examining the resulting stress field near the crack tip and the stress intensity 
factors. The results are given in Fig. 5. As the elements increase from eight to 20, the Kj (i = I, 
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II, III) alter only by less than 4%, 1%, and 1%, respectively. Thus, the use of a proper number of 
the singular element around the crack tip ensures the solution accuracy and convergency, if n > 8. 

5.3 Effect of Lamina Thickness and Delamination Location 

An important issue investigated in this study is the effect of laminate thickness, since 
transverse deformation and failure are affected significantly by lateral/transverse constraints in a 
composite laminate. In this paper, a delaminated [45/-45] s graphite-epoxy laminate with various 

b/h's are considered with b being kept constant. For an edge-delamination case, as shown in Fig. 

6, the results indicate that, for a composite with given laminate configuration, the change in ply 
thickness b/h does not alter the values of Kj significantly. The high value of Kjji relative to Kj 

indicate that mode-III crack growth and fracture may be the dominant failure mode. However, in 
a case of a delamination with short length, interaction of the free edge with the delamination 
becomes crucial. In the case of a center delamination in the composite laminate as shown in Fig. 7, 

it is found that the delamination crack growth becomes more dominated by mode I deformation as 
the ply thickness increases. The high value of Kj introduced by the increase in laminate thickness 

provides a clear basis for evaluation of failure mechanisms in thick-section composite laminates. 

5.4 Effect of Crack Length 

Fracture of a composite laminate resulting from delamination crack growth, is often 
observed in the experiment. Failure processes may experience stable, unstable crack growth or a 
combination of these. Thus, the effect of crack length on interlaminar crack-tip stresses may 

provide a rational basis of understanding the delamination growth behavior. The rapid increase in 
Kj and Gj for the center delamination cases in a graphite-epoxy composite is shown in Figs. 8-14. 

In contrast to the uniform stretching case, fracture of a delaminated composite subjected to bending 
may be catastropic, as all of the Kj and Gj increases monotonically as the value of a/b increases 

rapidly. 
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5.5 Effect of Laminate Layup 

Figures 15-21 illustrate the significant effect of fiber orientation on the crack-tip 
interlaminar stresses in [0/0] s graphite/epoxy thick composite laminates with various delammation 
lengths. For a center delamination case, the value of Kj varies appreciably with as depicted in 
Fig. 15. The Kj value reaches a plateau as 0 approaches 45°. As the crack length increases, a 
larger Kj is observed. A similar trend is shown in Fig. 16 for Kjj. Opposite tendency of the crack 
length effect is noticed. The important fracture mechanics parameter, Kjjj versus 0, is shown in 
Fig. 17. The value of Km increases first and then decreases with the ply angle 0. Note that 0 = 
37.5° is observed to introduce the maximum Kjjj for all of the crack lengths studied. For an edge 
delamination case, a similar behavior to the center delamination case is observed for the Kjij. 
However, the longer the delamination, the larger the Km. A similar behavior is given in Figs. 18- 
21 in which Gp Gjj and Gm as a function of ply angle 0 are shown. 
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6. Conclusions 

In this study, both a crack-tip asymptotic analysis and advanced numerical methods have 
been introduced to solve the delamination crack problem in a fiber composite laminate subjected to 
out-of-plane bending. The structure of the near-field solution has been determined by an 
eigenfunction expansion method. The full-field solution is determined by a conformal-mapping 
singular finite element analysis, based on the mathematical structure of the elasticity solution, 
numerical results are obtained for various cases with different lamination and crack parameters in 
the composite bending problem. Based on the methods developed and the results obtained, the 
following conclusions may be reached: 

1. The crack-tip asymptotic solution for a delaminated composite laminate under 
bending can be determined by an eigenfunction expansion method. 

2. A near-field singular element can be constructed on the basis of the asymptotic 
solution obtained and the conformal mapping method introduced by Stem [13]. 

3. In conjunction with the near-field singular element, a nonsingular composite 
bending element is also constructed based on the composite laminate elasticity 
solution for bending. 

4. The numerical methods developed lead to both the local crack-tip solution and 
and the far-field deformation and stress in a delaminated composite under out-of-plane 
bending. 

5. The effect of composite laminar thickness is found to be very significant on the 
delamination crack-tip stress intensities and associated energy release rates for a 
composite containing an interior delamination under external bending. However, much 
less appreciable influence is observed for a delamination located at the edge of the 
composite under nominal bending. 

6. The increase in laminar thickness may alter the failure mode of an interior delamination 
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under bending. Failure modes change from shear modes to the opening mode in a 
composite laminate subjected to out-of-plane bending. However, no appreciable change in 

the failure mode is found for an edge delamination case. 

7. In an angle-ply [0/ -0] s composite laminate containing a delamination in the 

upper interface, the fiber orientation 0 has very significant effects on the delami- 
nation crack-tip stress intensities and energy release rates. 

8. Significant delamination surface closure is found in composite specimens with 
both an edge delamination and a center interior delamination at the upper interface in the 
composite laminate under bending for the cases of 0 < 7.5°. Shear failure, thus, dominates 

the composite delamination crack growth behavior. 
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